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ABSTRACT 

Turbulence driven by magnetorotational instability (MRI) crucially affects the evolution of solid bodies in 
protoplanetary disks. On the other hand, small dust particles stabilize MRI by capturing ionized gas particles 
needed for the coupling of the gas and magnetic fields. To provide an empirical basis for modeling the co- 
evolution of dust and MRI, we perform three-dimensional, ohmic -resistive MHD simulations of a vertically 
stratified shearing box with an MRI-inactive "dead zone" of various sizes and with a net vertical magnetic flux 
of various strengths. We find that the vertical structure of turbulence is well characterized by the vertical mag- 
netic flux and three critical heights derived from the linear analysis of MRI in a stratified disk. In particular, 
the turbulent structure depends on the resistivity profile only through the critical heights and is insensitive to 
the details of the resistivity profile. We discover scaling relations between the amplitudes of various turbulent 
quantities (velocity dispersion, density fluctuation, vertical diffusion coefficient, and outflow mass flux) and 
vertically integrated accretion stresses. We also obtain empirical formulae for the integrated accretion stresses 
as a function of the vertical magnetic flux and the critical heights. These empirical relations allow to predict 
the vertical turbulent structure of a protoplanetary disk for a given strength of the magnetic flux and a given 
resistivity profile. 

Subject headings: dust, extinction — planets and satellites: formation — protoplanetary disks 



1. INTRODUCTION 

Planets are believed to form in protoplanetary gas 
disks. The standard scenario for planet formation con- 
sists of the following steps. Initially, submicron-sized dust 
grains grow into kilometer-sized planetesim als by coll i sional 
sticking and/or g r avitati o nal instability (ISafronoyl 119691 : 
iGoldreich & Wardf fl973t IWeidenschilling & Cuzzil [19931 
Planetesimals undergo further growth toward Moon-sized 
protoplanets thro ugh mutual collision assis ted by gravita- 
tional interaction ( Weth erill & Stewartlll989l) . Accretion of 
the dis k gas onto protoplanets leads to the formation of gas 
giants lzun o| [T980l iPollack et al] [1991 . while terrestrial 
planets form through the giant impacts of protoplanets after 
the gas disk disperse s by viscous accretion onto t he central 
star and other effects dChambers & Wetheriill998l) . 

Turbulence in protoplanetary disks plays a decisive role 
on planet formation as well as on disk dispersal. The 
impact of turbulence is particularly strong on the forma- 
tion of planetesimals since the frictional coupling of gas 
and dust particles governs the process. Classically, plan- 
etesimal formation has been attributed to the collapse of 
a dust sedimentary laye r by self-gravity dSafronovl Il969t 
IGoldreich & Ward 11973b and/or the collisio nal growth of 
dust grains ( IWeidenschilling & CuzzH [19931) . The pres- 
ence of strong turbulence is preferable for dust growth 
when the du st particles is so small that Coulomb repulsion 
is effective dOkuzum 1 120091 lOkuzumi et aT] l20llh . How- 
ever, strong turbulence acts against the growth of macro- 
scop ic dust aggregates sinc e it makes their collis ion disrup- 
tive (|\^idenschnnngl [1984]; [Johansen]et]n][2008|). Further- 
more, turbulence causes the diffusion of a dust sedimen- 
tary layer, making planetes imal formation via gra vitational 
instability difficult as well (IWeidenschi lling 1984). Turbu- 



lence is also known to concentrate dust particles of particular 
sizes, but its relevance to planetesimal f ormation via grav- 
itational instability is still under debate dCuzzi et al.l 1200 11 
l200ll20loHPan et alj|201lb . More recently, it has been sug- 
gested that two-fluid instability of gas and dust can produce 
dust clumps with density high enough for gravitational col- 
lapse, but successful dust coagulation to macroscopic sizes 
seems to be still required for this mechanism to become 
viable (|Youdin & Goodma n 2005; Johansen & Youdi"n1l2007l : 
Uohansen et alJl2007tlBai & Stonell2010l) . Besides, turbulence 
also affects planetesimal growth as turbulent density fluctua- 
tions gravitationally interact with planetesimals and can raise 
their random velocities above the escape velocity dlda et alj 
120081: iNelson & Gressej||2010l) . The fluctuating gravitational 
field can even cause random orbital migration of p rotoplanets 
dLaughlin et al.H2004 INelson & Papaloizoull2004l) . Thus, to 
understand the growth of solid bodies in various stages, it is 
essential to know the strength and spatial distribution of disk 
turbulence. 

Interestingly, the evolution of solid bodies is not only af- 
fected but also affects disk turbulence. The most viable mech- 
anism for generat ing disk turbulence is the magnetorotational 
instability (MRI; iBalbus & Hawley||1991l) . This instability 
has its origin in the interaction between the gas disk and mag- 
netic fields, and therefore requires a sufficiently high ioniza- 
tion degree to operate. Importantly, whether the MRI oper- 
ates or not in each location of the disk is strongly depen- 
dent on the amount of small dust grains because they ef- 
ficiently capture i onized gas particles and thus reduce the 
ionization degr ee dSano et aUboOOt lllgner & Nelson! 120061; 
Okuzumi 2009). This implies that dust and MRI-driven tur- 
bulence affect each other and thus evolve simultaneously. 

The purpose of this study is to present an empirical ba- 
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sis for studying the coevolution of solid particles and MRI- 
driven turbulence. It is computationally intensive to simulate 
the evolution of dust and MRI-driven turbulence simultane- 
ously, since the evolutionary timescale of solid bodies is gen- 
erally much longer than the dynamical timescale of the tur- 
bulence. For example, turbulent eddi es grow and decay on a 
times cale of one orbital period (e.g.. iFromang & Papaloizoul 
120061) . while dust particles grow to macroscopic sizes and 
settle to the midplane spen d ing 100-1000 orbital p eriods 
(e.g.. iNakagawa et al.l 119811: iDullemond & Dominikl 12005b 
iBrauer et al.l 120081) . However, this also means that MRI- 
driven turbulence can be regarded as quasi-steady in each evo- 
lutionary stage of dust evolution. Motivated by this fact, we 
restrict ourselves to time-independent ohmic resistivity, but 
instead focus on how the quasi-steady structure of turbulence 
depends on the vertical profile of the resistivity. 

To characterize the vertical structure of MRI-driven tur- 
bulence, we perform a number of three-dimensional MHD 
simulations of local stratified disks including resistivity and 
nonzero net vertical magnetic flux. Inclusion of a nonzero 
net vertical flux is im portant as it determines the satura- 
tion level of turbulence ( | Hawlev et al.ll 1 995b ISano et al.l[2004t 
ISuzuki & Inutsukal I2009L see also our Section SJ. Simi- 
lar simulatio ns have been done i n a number of previous 
studies (e.g.jMiller & Stonell2000t ISuzuki & Inutsukall2009t 



Oishi & Mac Low! 120091: ISuzuki et al.l l2010t iTurneretal 
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201 11) . One important difference between our study and pre- 
vious ones is that we focus on general dependence of the sat- 
urated turbulent state on the model parameters such as the re- 
sistivity and net magnetic vertical flux. 

Our modeling of MRI-driven turbulence follows two steps. 
In the first step, we seek scaling laws giving the rela- 
tions among turbulent quantities. We express the relations 
as a function of the vertically integrated accretion stress, 
which is the quantity that determines the rate at which 
turbulent energy is extracte d from the differential rotation 
dBalbus & PapaloizoulfT999l) . As we will see, excellent scal- 
ing relations are obtained if we divide the integrated stress 
into two components that characterize the contributions from 
different regions in the stratified disk (which we will call the 
"disk core" and "atmosphere'') In the second step, we find 
out empirical formulae that predict the vertically integrated 
stresses as a function of the resistivity profile and vertical 
magnetic flux. 

The plan of this paper is as follows. In Section [2] we de- 
scribe the method and setup used in our MHD simulations. 
In Section|3] we introduce "critical heights" derived from the 
linear analysis of MRI in stratified disks. As we will see later, 
these critical heights are useful to characterize the turbulent 
structure observed in our simulations. We present our sim- 
ulation results in Section @] and obtain scaling relations and 
predictor functions for the quasi-steady state of turbulence in 
Section In Section [6] we simulate dust settling in a dead 
zone to model the diffusion coefficient for small particles as 
a function of height. Effects of numerical resolutions on our 
simulation results are discussed in Section|7] Our findings are 
summarized in Section[8] 

2. SIMULATION SETUP AND METHOD 

In this section, we describe the setup and method adopted 
in our stratified resistive MHD simulations. 

2.1. Setup 



O ur MHD simulations adopt the shearing box approxima- 
tion (lHawlev et al.lll995l) . We consider a small patch of disk 
centered on the midplane of an arbitrary distance from the 
central star, and model it as a stratified shearing box corotat- 
ing with the angular speed O at the domain center. We use the 
Cartesian coordinate system (x,y,z), where x, y, and z stand 
for the radial, azimuthal, and vertical distance from the do- 
main center. In addition, we assume that the gas is isothermal 
throughout the box; thus, the sound velocity c s of the gas is 
constant in both time and space. 

2.1.1. Initial Conditions 

For the initial condition, we assume that the gas disk is ini- 
tially in hydrostatic equilibrium and is threaded by uniform 
vertical magnetic field B z q. The assumption of the hydrostatic 
equilibrium leas to the initial gas density profile 



p= poexp 



2h 2 J 



where po is the initial gas density at the midplane and 



h = 



(1) 



(2) 



is the pressure scale height[] The ratio between the initial 
midplane gas pressure Po = A) c .? ar, d me initial magnetic pres- 
sure B 2 /87T defines the initial plasma beta 



fa 



&np c 2 
B 2 



(3) 



In this paper, the strength of the initial magnetic flux will be 
referred by /3~q rather than B z q. 

2.1.2. Resistivity Profile 

The main purpose of this study is to see how the turbulence 
depends on the vertical profile of the ohmic resistivity. We 
adopt a simple analytic resistivity profile based on the follow- 
ing consideration. For fixed temperature, the resistivity is in- 
verse ly proportional to the ionization degree ( Blaes & Balbusl 
119941) . In protoplanetary disks, the ionization degree at each 
location is determined by the balance of ionization (by, e.g., 
cosmic rays and X-rays) and recombination (in the gas phase 
and on grain surfaces). Detailed structure of the resistivity 
profile depends on what processes dominate the ionization 
and recombination. However, a general tendency is that the 
ionization degree decreases toward the midplane of the disk, 
because the ionization rate is lower as the column depth is 
greater and because the reco mbination rate is higher as the gas 
density is higher (see, e.g., ISano et ail|2000l) . Based on this 
fact, we give the resistivity profile rj{z) such that 77 increases 
as z decreases. To be more specific, we adopt the following 
resistivity profile 



?? = ?7midexp 



2h% 



(4) 



where rj m a is the resistivity at the midplane and h n is the scale 
height of 77. 

Equation (01 satisfies the important property of realistic re- 
sistivity profiles mentioned above. Furthermore, as shown in 

1 Note that the "gas scale height" is often defined as H = \/2h = \/2c s /Q 
in the literature on stratified MHD simulations. 
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Table 1 

Model Parameters and Initial Critical Heights 
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Appendix, Equation (01 exactly reproduces the vertical resis- 
tivity profile of a disk in some limited cases. However, it will 
be useful to examine possible influences of limiting rj to Equa- 
tion dH) . In order to do that, we a lso consider a resistive profile 
used in lFleming & Stond (120031) . 



?7fso3 = ?7oexp 



z 



exp 




(5) 

which is characte rized by two parameters 770 and So /Scr (see 
Equation (10) of IFleming & Stonell2003l) . Physically, Equa- 
tion (0 corresponds to the resistivity profile when the ioniza- 
tion degree is determined by the balance between cosmic -ray 
ionization and gas-phase recombination (see also Appendix). 

We construct 17 simulation models using Equations and 
(O. 16 models are constructed from Equation © with var- 
ious sets of the parameters (77^, h v , (3 z q). Table Q] lists the 
parameters adopted in each run. The model "Ideal" assumes 
zero resistivity throughout the simulation box. Models X0- 
X3 are defined by the same value of but different values 
of h v . The difference between X, Y, and W models is in the 
value of 77 m j(j. The initial plasma beta is taken to be 3 x 10 5 
in all models except Xla-Xld. In addition, we construct one 
model from Equation (O with r) m iA = 7yoexp(So/4EcR) = 0.01 
and Eo/Scr = 34.7. This set of parameters correspond s to 
the "larger dead zone" model of IFleming & Stoiiel (12003b de- 
fined by ReM,mid = 100 and Re M . =2v ^/, = 5.6 x 10 6 , where 
Re« = c s h/rj is the magnetic Reynolds number with the typi- 
cal length and velocity set to be h and c s , respectively. We re- 
fer to the model with these parameters as model FS03L. The 
initial plasma beta in the FS03L model is taken to be 3 x 10 s . 
Note that our FS 03L run is not exactly th e same as the larger 
dead zone run o f IFleming & Stoiiel ( 120031) since they assumed 
zero net vertical magnetic flux. 

2.2. Method 

We solve the equat ions of the isothermal re sistive MHD us- 
ing the ZEUS code dStone & Normanll 19921) . The domain is 
a box of size 2^/2h x 8\/2h x 10\/2h along the radial, az- 
imuthal, and vertical directions, divided into 40 x 80 x 200 



grid cells. For all runs, the radial and azimuthal boundary con- 
ditions are taken to be shearing-periodic and periodic, respec- 
tively. Therefore, the vertical magnetic flux is a conserved 
quantity of our simulations. The vertical boundary condition 
for all runs except Xld is the standard ZEUS outflow condi- 
tion, where the fields on the boundaries are computed from 
the extrapolated electromotive forces; only for run Xld, we 
assume for numerical stability that the magnetic fie lds are ver- 
tical o n the top and bottom boundaries as done bv lFlaig et al.l 
(120101) . We have checked that the results hardly depend on the 
choice of the two types of boundary conditions. Note that the 
radial and azimuthal components of the mean magnetic fields 
are not conserved due to the outflow boundary condition. We 
also added a small a rtificial resistivity n ear the boundaries for 
numerical stability dHirose et al.l 12009). A density floor of 
10~ 5 po is applied to prevent high Alfven speeds from halting 
the calculations. 

3. CHARACTERISTIC WAVELENGTHS AND CRITICAL 
HEIGHTS 

As we will see in the following sections, it is useful to 
analyze the simulation results using the knowledge obtained 
from the linear analysis of MRI. In this subsection, we intro- 
duce several quantities that characterize the linear evolution 
of MRI. 

3.1. Characteristic Wavelengths of MRI 

According to t he local linear analysis of MRI including 
ohmic resistivity dSano & Mivamal[l999l) . the wavelength of 
the most unstable MRI mode can be approximately expressed 
as 

Alocal ~ max{ Aideal , Ares } , (6) 



where 



and 



Aideal = 2-7T 



VAz 
O 



VAz 



(7) 
(8) 



are the characteristic wavelengths of MRI modes in the ideal 
and resistive MHD limits, respectively, and va z = B z j y/AiFp is 
the vertical component of the Alfven velocity. Equation (O 
can be written as Ai oca i ~ Aid ea imax{l, A" 1 }, where A is the 
Elsasser number defined by 



Az 
7]Q' 



(9) 



The Elsasser number determines the growth rate of the MRI. 
If A > 1, ohmic diffusion does not affect the most unstable 
mode lying at A = Aid ea i> an d the local instability occurs rapidly 
at a wavelength Ai oca i w Aideal and at a rate sa fl. If A < 1, 
ohmic diffusion stabilizes the most unstable mode, and the 
local instability occurs at a longer wavelength Ai oca i w A res = 
A -1 Aideal and at a slower rate ps A _1 fl We will refer to the 
former case as the "ideal MRI," and to the latter case as the 
"resistive MRI." 

Figure Q] schematically illustrates how Aideal and A res vary 
with height \z\. In general, Aideal grows toward higher \z\ be- 
cause the Alfven speed vaz increases as the density decreases. 
By contrast, A res grows toward lower |z| because A res is in- 
versely proportional to va z and because rj incre ases with de- 
creasing \z\ (see the discussion in Section l2.1.2l ). 
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Figure 1. Schematic illustration showing the vertical structure of a stratified 
protoplanetary disk with vertical magnetic fields. The horizontal axis shows 
the distance z from the midplane, while the vertical axis shows the charac- 
teristic wavelengths A[d ea i (solid curve) and A res (dashed curve) of MRI at 
each z as well as the gas scale height h (dotted line). The set of three ra- 
tios AyealA, A r c S //2, and Aidcal/Ares = A defines four layers. At \z\ > /lideal 
(A > 1 and Ajd ca i > h), MRI is stabilized due to the weak gas pressure com- 
pared to the magnetic tension. At h\ < \z,\ < ftideal (A > 1 and A[d ea i < 
dark gray region), MRI operates without affected by ohmic dissipation. At 
/'res < |z| < (A < 1 and A res < h; light gray region), ohmic dissipation is 
effective and MRI operates only weakly. At \z\ < Arcs (A < 1 and Arcs > h), 
ohmic dissipation perfectly stabilizes MRI. We refer to the regions \z\ < ^ideal 
and \z\ > /lidcal as me "disk core" and "atmosphere," respectively. 

The global instability of a stratified disk can be described 
in term s of the local analysis. As shown by ISano & Mivamal 
( 119991) . the gas motion at height z is unstable if the local un- 
stable wavelength Ai oca i is shorter than the scale height of the 
disk, i.e., 



max{A idea i,A res } < h. 



(10) 



3.2. Critical Heights 

With the global instability criterion (Equation ( flOl i) together 
with the vertical dependence of Ajdeai and A res , we can define 
three different critical heights for a stratified disk. 



1 . The first one is /ijdeai defined by 

Aidealfe = ^ideal) = h, 



(ID 



or equivalently, f3„(z = /ijdeai) = oV 2 , where f3 z (z) = 
&irp(z)c 2 s /B 2 z (z). At \z\ > /iideai Wz £ 8tt 2 ), MRI does not 
operate because the wavelengt hs of the unstable mode s 
exceed the disk thickness ~ h (ISano & Miyamairi999l) . 
We refer to the region \z\ ^ /Jideai as the "atmosphere" 
and to the region \z\ ^ /z,deai as the "disk core." 



2. The second one is defined by 

A(z = A A ) = l, 



(12) 




-4-2 2 4 
height z/h 

Figure 2. Vertical profiles of the ohmic resistivity r) (upper panel) and the 
initial Elsasser number A (= o (lower panel) for models X0 (blue dot-dashed 
curve), XI (blue dotted curve), X2 (blue dashed curve), X3 (blue solid curve), 
Y3 (red curve), W3 (green curve), and FS03L (black curve). 

regard as zero when a dead zone is absent. This is 
the case for model Ideal. 



3. The third one is /z les defined by 

Ares(z=^res)=/J- 



(13) 



Ohmic diffusion allows the resistive MRI to operate at 
/ 1 res ;$ \z\ < /ja- At \z\ < li,es, ohmic diffusion stabilizes 
all the unsta ble MRI modes. Note that some previous 
studies (e.g.. lGammidll996tlSano et al.ll2000l) used the 
terminology "dead zone" for the region \z\ ^ h K s rather 
than \z\ ^ Iia- In fact, as we will see in Section|4j the 
set of h\ and h KS best characterizes our dead zone. We 
regard /z les as zero when A res is less than h at all heights. 
This is the case for Y models. 

The critical heights in the initial state (Mdeaio> / 1 A.o> an d 
^res.o) are shown in Table 1 for all of our 17 simulations. Us- 
ing Equations (Q]) and one can analytically calculate the 
initial critical heights for models except FL03L as 



h 



ideal, ' 



hA,0 = 



2 In 



21nAo' 

i+(W 2 



1/2 



1/2 



'res,0 : 



or equivalently, Ai dea i(z = h a) = A les (z = h^)- The layer 
h\ < \z\ ^ ^ideai is the so-called active layer, where where 
MRI operates without affected by ohmic diffusion nor 
gas stratification. The region \z\ ^ h\ is what we call 
the "dead zone," where ohmic diffusion stabilizes the 
most unstable ideal MRI mode. For convenience, we 



21n(87r 2 /3- 1 A^ 2 
\+2(h/h v ) 2 

A f) = 2 ^ 



h, 



^midAo 

is the initial Elsasser number at the midplane. 



(14) 



(15) 



(16) 



(17) 
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Figure 3. Horizontally averaged Maxwell stress wu (upper panel) and 
density-weighted velocity dispersion p<5v 2 (bottom panel) as a function of 
time ( and height z for model X 1 . The solid, dotted, and dashed lines show 
the critical heights z = /Jideal* ' ! A> and ''res, respectively. 

Figure |2] shows the vertical profiles of the resistivity 77 and 
the initial Elsasser number A f= o for some of our models. The 
initial midplane Elsasser number Ao and initial critical heights 
(/zideai,o,^A,o,^res.o) are listed in Table[TJfor all models. As one 
can see from Table [TJ and the lower panel of Figure [2] models 
labeled by the same number are arranged so that they have 
similar values of h\fi. 

For turbulent states, we evaluate va z m the Elsasser number 
and the characteristic wavelengths as (B 2 j\Tfp) x l 2 , where the 
overbars denote the horizontal averages. 

4. SIMULATION RESULTS 
4. 1 . The Fiducial Model 

We select model XI as the fiducial model to describe in 
detail. Figure [3] shows how MRI-driven turbulence reaches 
a quasi-steady state in run XI. The upper and lower pan- 
els plot the horizontal averages of the Maxwell stress Wm = 
—SB x SB y /A-K and the density-weighted velocity dispersion 

~pSv 2 , respectively, as a function of time t and height z. The 
solid, dotted, and dashed lines are the loci of the critical 
heights Ziideai, h\, and /i res , respectively. As seen in the fig- 
ure, a quasi-steady state is reached within the first 40 orbits. 
The critical height //ideal measured in the quasi-steady state is 
slightly lower than that in the initial nonturbulent state. This is 
because the ideal MRI wavelength Aideai oc va z is increased by 
the fluctuation in the vertical magnetic field, SB 7 :. By contrast, 
h\ and /z res are almost unchanged, because the fluctuation of 
the magnetic field is suppressed in the dead zone. 

Figure |4] shows the vertical structure of the disk averaged 
over a time interval 175 < f2?/(27r) < 350. The dark and light 
gray bars in each panel indicate the heights where ideal and 
resistive MRIs operate, respectively (see also Figure [TJ. The 
brackets (• ••) denote the averages over time and horizontal 
directions. 

In Figure Ufa), we compare the averaged gas density (p) 
with the initial density given by Equation (UJ). We see that 
the density is almost unchanged in the disk core (|z| < /zideai) 
but is considerably increased in the atmosphere (\z\ > Aideai)- 



This is because the magnetic pressure is negligibly small in 
the disk core but dominates over the gas pressure in the atmo- 
sphere. Figure |4j a) also shows the amplitude of the density 
fluctuation, (Sp 2 ) 1 ! 2 . As one can see, the density fluctuation 
is small ((Sp 2 ) 1 ! 2 < (p)) except at |z| > /i idea ]. 

The magnetic activity in the disk can be seen in Fig- 
ure |4|b), where the vertical profiles of the magnetic energies 
((SB 2 )/4tt, (B 2 ,)/4tt, and (B 2 )/4tt) and Maxwell stress (w M ) 
are plotted. One can see that these quantities peak near the 
outer boundaries of the active layers, \z\ ~ hide-A- This is be- 
cause the largest ch annel flows develop at lo cations where 
Aideai ~ h (see, e.g., ISuzuki & InutsukaTl2009t) . In the dead 
zone, ohmic dissipation suppresses the fluctuation in the mag- 
netic fields, (SB 2 ), leaving the initial vertical field ((B 2 ) ps B 2 q ) 
and coherent toroidal fields ((B 2 ) sa (By) 2 ) generated by the 

differential rotationQ 

Figure Hfc) shows the density-weighted velocity disper- 
sions (p)(8v 2 ) and (p)(Sv 2 ) and the Reynolds stress (wr) = 
(pSv x Svy), These quantities characterize the kinetic energy in 
the random motion of the gasQ Comparing Figures |4|b) and 
(c), we find that the drop in these quantities in the disk core 
is not as significant as the drop in (SB 2 ) and (%). This is 
an indication that sound waves gene rated in the active layer s 
penetrate deep inside the dead zone (iFleming & Stondl2003l) . 
Furthermore, we find that (p)(Sv 2 ) is approximately constant, 
i.e., the velocity dispersion (i5v 2 ) is inversely proportional to 
the mean density (p), in the disk core. This means that the 
kinetic energy density of fluctuation is nearly constant in the 
disk core. This is another indication of sound waves, because 
the amplitude of the velocity fluctuation Sv is generally pro- 
portional to 1/ Jp for freely propagating sound waves. The 

root-mean-squared random velocity (Sv 2 ) 1 ! 2 j s shown in Fig- 
ure|4jd). The random velocity is subsonic in the disk core and 
exceeds the sound speed only in the atmosphere. 

An indication of freely sound waves can be also found in 
the density fluctuation. Shown in Figure Hfe) is the mean 
squared density fluctuation (Sp 2 ) divided by the mean den- 
sity (p). Since (Sp 2 ) 1 ! 2 < (p), the quantity (Sp 2 )/(p) ps 
(Sp 2 /p) is approximately proportional to the thermal energy 
density of fluctuation, (c 2 Sp 2 /2p). In the disk core, we see 
that (Sp 2 ) I (p) is roughly constant along the vertical direc- 
tion, meaning that the amplitude of the density fluctuation, 
(Sp 2 ) 1 ' 2 , is proportional to the square root of the mean den- 
sity (p). The similarity between (p)(Sv 2 ) and (8p 2 )/(p) is 
peculiar to sound waves, for which 8v/c s ~ Sp/ p. 

Figure |4jf) displays the profile of the vertical mass flux 
(pv z ). In the atmosphere, the vertical flux is outward, i.e., 
(pv z ) > at z > Aideai and (pv z ) < at z < -^ideai- This outflow 
results from the breakup of large channe l flows at the outer 
boundaries of the active la yers, |z| w /ijdeai d Suzuki & Inutsukal 
120091: ISuzuki et al.l 120101) . The vertical mass flux reaches 
a constant value at |z| > 5h. This fact allows us to mea- 
sure awell-defined outflow flux for each simulation (see Sec- 
tion [5T3l . 

4.2. Model Comparison 

2 In our simulations, ohmic resistivity is not high enough to remove shear- 
generated, coherent toroidal fields. In this sense, o ur dead zone is an "undead 
zone" in the terminology of Turner & Sano 1 2008). 

3 In the disk core (|z| < /lideal), (p) {<5v 2 } is approximately equal to (p5v 2 ) 
since the density fluctuation is small (see Figure|4ja)). 
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Model XI Model XI Model XI 




height z/h height z/h height z/h 

Figure 4. Vertical profiles of various quantities averaged over x—y planes and over a time interval 175 orbits < t < 350 orbits for model XI. (a) Initial (dashed 
curve) and time-averaged (solid curve) gas densities, and amplitude of the density fluctuation (dotted curve), (b) Magnetic energies {SB 2 )/4n (solid curve), 
(By) /4tt (dot-dashed curve), and (B 2 ) /Air (dashed curve), and Maxwell stress (wm) (dotted curve) normalized by the initial midplane gas pressure Po = Pqc 2 . (c) 
Density-weighted velocity dispersions (p) (<5v 2 ) (thick solid curve) and (p) (<5v? ) (dashed curve), and Reynolds stress (wr) (dotted curve), (d) Root-mean-squared 
random velocity (i5v 2 ) 1 I 2 . (e) (Sp 2 ) / (p) , which represents the thermal energy of fluctuation, (f) Vertical mass flux (pv : ) , plotted for positive (solid curve) and 
negative (dashed curve) values. The regions shaded in dark and light gray indicate where ideal and resistive MRIs operate, respectively (see also Figure[TJ. 




Figure 5. Vertical profiles of temporally and horizontally averaged Maxwell stress (wm) (red curves) and density-weighted velocity dispersion (p)(&v 2 ) (black 
curves) normalized by the initial midplane gas pressure Po = poc 2 for /3 z q = 3 X 10 5 models. The dark gray bars indicate where MRI operates without affected by 
ohmic resistivity (/ia < |z| < /iideal), while the light gray bars show where MRI operates but is weakened by ohmic resistivity (/i res < \z\ < h^). 
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We now investigate how the vertical structure of turbulence 
depends on the resistivity profile and vertical magnetic flux. 

Figure |5] displays the temporal and horizontal averages of 
the Maxwell stress (wm) and the density-weighted velocity 
dispersion (p) (8v 2 ) as a function of z for various /3-o = 3 x 10 5 
models. As in Figures[T]and|4] the dark and light gray bars in 
each panel indicate the heights where ideal and resistive MRIs 
operate, respectively. 

The effect of changing the size of the dead zone can be seen 
in Figure |3 a), where models XI, X2, and X3 are compared. 
These models are characterized by the same values of /?-o and 
Ao but different values of h n . For all the models, (wm) sharply 
falls at z ~ /ires, meaning that /i res well predicts where the resis- 
tivity shuts off the magnetic activity. By contrast, (p) (Sv 2 ) ex- 
hibits a flat profile at \z\ < /jjdeai with no distinct change across 
z = h KS nor z = ht\. The only clear difference is the value of 
(p) (8v 2 ) in the disk core, i.e., the value is lower when the dead 
zone is wider. Note that (p) (Sv 2 ) decreases more slowly than 
the column density of the active layers (h\ < \z\ < /Jideai)- For 
example, the active column density in model XI is 20 times 
smaller than that in model X3. However, the midplane value 
of (p) (5v 2 ) in the former is only five times smaller than that 
in the latter. This suggests that even a very thin active layer 
can provide a large velocity dispersion near the midplane. 

Interestingly, the vertical structure of turbulence depends on 
the critical heights (/ijdeai, h\, and h les ) but are very insensitive 
to the details of the resistivity profile. This can be seen in Fig- 
ure HJb), where we compare runs with similar critical heights 
(runs X3, W3, and FS03L). We see that these models pro- 
duce very similar vertical profiles of (wm) and (p) (Sv 2 ) even 
though they assume quite different resistivity profiles (see the 
upper panel of Figure[2]i. This suggests that the vertical struc- 
ture of turbulence is determined by the values of the critical 
heights. 

Importance of distinguishing hj\ and h KS is illustrated in 
Figures|51c) and (d). These panels compare five models (Yl- 
4 and Ideal) in which the resistive MRI is active at the mid- 
plane, i.e., /z res = 0. Figure|3c) shows models with Iia > Q.5h. 
We see that the profiles of (p) (Sv 2 ) and (wm) are very similar 
for the three models. This implies that the vertical structure is 
determined by the value of /i res when h\ > 0.5h. Figure |5Jd) 
shows what happens when h\ falls below 0.5h. Model Ideal is 
clearly different from the other models. Model Y4 (Iia = 0.4) 
is interesting because it exhibits both features. In the lower 
half of the disk (z < 0), the Maxwell stress behaves as in the 
other Y models. In the upper half (z > 0), however, the profile 
of (wm) is closer to that in model Ideal. 

Next, we see how the saturation level of turbulence depends 
on the vertical magnetic flux. Figure|6]shows the vertical pro- 
files of (wm) and (p) (Sv 2 ) for three runs with different values 
of /3,o (Xlb, XI, and Xld). We see that these values increase 
with decreasing /3-q. The peak value of (wm) is approximately 
10" 4 P , 10 _3 P , and lO" 2 ^ for runs Xlb, XI, and Xld, re- 
spectively (Pq = Poc 2 is the initial midplane gas pressure). This 
indicates a linear scaling between the turbulent stress and /3~o ■ 

5. SCALING RELATIONS AND PREDICTOR 
FUNCTIONS 

Now we seek how the amplitudes of turbulent quantities de- 
pend on the vertical magnetic flux and the resistivity profile. 
We do this in two steps. First, we derive relations between the 
amplitudes of turbulent quantities and the vertically integrated 
turbulent stress. We then obtain empirical formulae that pre- 
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Figure 6. Vertical profiles of temporally and horizontally averaged Maxwell 
stress (wm ) (red curves) and density- weighted velocity dispersion (p) (<5v 2 } 
(black curves) for models with different P z q. The dotted, solid, and dashed 
curves correspond to models Xlb, XI, and Xld, respectively. 

diet the integrated stress as a function of the vertical magnetic 
flux and the resistivity profile. 

5.1. Scaling Relations between Turbulent Quantities and 
Vertically Integrated Accretion Stresses 

The ultimate source of the energy of turbulence is the 
shear motion of the background flow. The accretion stress 
Wxy = wr + wm determines the rate at which the free energy 
is extracted. Therefore, we expect that the accretion stress is 
related to the amplitudes of turbulent quantities, such as the 
gas velocity dispersion and outflow mass flux. 

To quantify the rate of the energy input in the simulation 
box, we introduce the effective a parameter 



= J(w„)dz 
~ £c 2 ' 



(18) 



where £ = J (p)dz is the gas surface density. In the classical, 
one-d imensional viscous disk theory dLynden-Bell & Pringld 
119741) . the parameter a is related to the turbulent viscosity 
t'trub as i/ tru b = (3/2)ac 2 /fi, where the prefactor 3/2 comes 
from the slope of the Keplerian rotation. Thus, a also charac- 
terizes the vertically integrated mass accretion rate. 

As we will see below, it is useful to decompose a as a = 
cw + aatm, where 



and 



Or 



aatm ! 



|z|<ftideal 



[w xy )dz 



(19) 



(20) 



are the contributions from the disk core (|z| < /ijdeai) and atmo- 
sphere (|z| > /iideai), respectively. Table|2]shows the values of 
ot, a core , and a a tm as well as the time-averaged critical heights 
{hide-dh liA, hres) for all our simulations. 

5.1.1. Velocity Dispersion 

Random motion of the gas crucially affects the growth 
of dust particles as it enhances the collision velocity be- 
tween the particles v ia friction forces iWeidenschilUngiri984t 
iJohansen et alj|2008l) . Here, we seek how the velocity disper- 
sion (Sv 2 ) is related to the integrated accretion stress. 

First, we focus on the velocity dispersion at the midplane, 
(^v 2 )mid. Figure|7ta) shows (Sv 2 ) m id versus the total accretion 
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Table 2 

Time-averaged Properties of MHD Simulations 



Model 


''ideal 

h 


"A 

h 


h 

_!}2: 
h 


ft 

10~ 3 


10~ 3 


a 

— 
10~ 3 


( f)\^\.~;,i 
\ u v / mid 

10~ 3 cf 


\ Z / nua 

(5v 2 ) m id 


\ U H i mid 


10" 5 (p)midCj 


Ideal 


2.4 


0.0 


0.0 


18 


12 


6.0 


15 
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4.7 


XO 


3.3 
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0.13 
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2.5 


XI 


3.2 


2.5 


2.2 


2.0 


0.42 


1.5 
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Figure 7. Gas velocity dispersion at the midplane, (<5v )nM, for all runs 
presented in this study. Panels (a) and (b) plot the data versus a and a CO re, 
respectively. The symbols correspond to models Ideal (circle), X0-X3 (open 
squares), Y1-Y4 (triangles), W1-W3 (crosses), Xla-Xld (filled squares), 
and FS03L (plus sign). The lines show the best linear fits (Equation )2U for 
panel (b)). 



stress a for all our runs. The value of (<5v 2 ) m jd for each run is 
listed in Table [2] One can see a rough linear correlation be- 
tween the velocity dispersion and the accretion stress (for ref- 
erence, a linear fit (<5v 2 ) m jd = 0.25ac 2 is shown by the dashed 
line). However, detailed inspection shows that (<5v 2 ) m id de- 
creases more rapidly than a as the dead zone increases in size. 
As found from Table 2, this is because the contribution from 
the atmosphere, a atm , is insensitive to the size of the dead 
zone in the disk core. In Figure |7{b), we replot the data by 
replacing a with the accretion stress in the disk core, a C0le - 
Comparison between Figures|7Ja) and (b) shows that (<5v 2 ) m jd 



more tightly correlates with a core rather than with a. We find 
that the data can be well fit by a simple linear relation 



(<5v 2 ) m id = 0.78a core c 2 , 



(21) 



which is shown by the solid line in Figure Hb). This result 
indicates that the accretion stress in the atmosphere does not 
contribute to the velocity fluctuation near the midplane. 

Once (5v 2 ) m id is known, it is also possible to reproduce the 
vertical profile of the velocity dispersion. For the disk core 
(\z\ < /iideaiX we already know that (i5v 2 ) is inversely propor- 
tional to the mean gas density (p) and that (p) hardly deviates 
from the initial Gaussian profile. From these facts, we can 
predict the vertical distribution of (<5v 2 ) as 



(6v 2 )^(6v 2 ) 



mid " 



(p) 



mid 



:0.78a r 



(P) 
e c 2 exp^ 



(£v 2 ) mid expf 



2h 2 J 



z 2 

2k 2 



(22) 



where Equation (l2Tb has been used in the final equality. In 
Figure [8] we compare the vertical profiles of the random ve- 
locity (8v 2 ) 1 ' 2 directly obtained from runs Ideal, XI, and Xla 
with the predictions from Equation ( 1221 ). where the values 
of a C ore are taken from Table 2. We see that Equation (1221 
successfully reproduces the vertical profiles of (5v 2 ) l ' 2 in the 
disk core. We remark that Equation ( 1221 greatly overestimates 
the velocity dispersion at |z| h-^i, where the gas density 
can no longer be approximated by the initial Gaussian profile 
(see FigureHfa)). 

5. 1 .2. Density Fluctuation 

Density fluctuations generated by MRI-driven turbulence 
gravitationally interact with planetesimals and larger solid 
bodies, affecting t heir collisional and orbital evolution in pro- 
toplanetarv disks (lLaughlin et a l. 2004; Nelson & Papaloizoul 



l200^lNeTson & Gresselll2010tlGressel et al.ll2011l) . Here, we 
examine how the amplitude of the density fluctuations is de- 
termined the ve rtically integrated accretion stress. 

As in Section 15.1.1 1 we begin with the analysis of the den- 
sity fluctuations at the midplane, (^p 2 )^- We find from Ta- 
ble|2]that (Sp 2 ) m a more tightly correlates with a core than with 
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Figure 8. Vertical profiles of the random velocity (<5i' 2 }'/ 2 directly obtained 
from MH D si mulations (black curves), compared with the predictions from 
Equation {22} with a C orc taken from Table[2](gray curves). The solid, dashed, 
and dotted curves correspond to runs Ideal, XI, and Xla, respe ctiv ely. The 
predicted profiles are plotted only at \z\ < Aideal, where Equation (22) is valid. 




Figure 9. Gas density fluctuation at the midplane, {<5p 2 ) m id/ (p)mid' versus 
"core f° r a " mns presented in this study. The symbols correspond to models 
Ideal (circle), X0-X3 (open squares), Y1-Y4 (triangles), W1-Y3 (crosses), 
Xla-Xld (filled sq uare s), and FS03L (plus sign). The line shows the best 
linear fit (Equation )23t ). 



Figure [9] shows (<5p 2 )mid/(p)mid versus ct core for all runs. 



The best linear fit is given by 



(5p 2 



mid 



■0A7a mie (p) 



2 

midi 



(23) 



which is shown by the solid line in Figure [9] If we use this 
equation with Equation d2lT) . we can also obtain the rela- 
tion between the velocity dispersion and density fluctuation, 
(<5v 2 ) mid /c 2 = 1.7((5p 2 ) mid / (p);^. This is consistent with the 
idea that the fluctuations near the midplane are creat ed b y 
sound waves, for which Sv/c s ~ Sp/p (see also Section l4Tt . 

As shown in Section |4~T1 (dp 2 ) is roughly proportional to 
(p) along the vertical direction in the disk core. Hence, if 

(^/° 2 )mid i s § iven > one can reconstruct the vertical profile of 
the density fluctuation in the disk core according to 



(5p 2 ) « (<5/9 2 )mid 



(P) 



(p)mid 

0.47a core (p) 2 nid exp 



z 2 



(5p 2 )midexp 

where we have used (p) ss (p) m idexp(-z 2 /2/i 2 ) and Equa 



(24) 
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Figure 10. Outflow mass flux m w normalized by (p) m idCj for all runs pre- 
sented in this study. Panels (a) and (b) plot the data versus a and « a tm, re- 
spectively. The symbols correspond to models Ideal (circle), X0-X3 (open 
squares), Y1-Y4 (triangles), W1-W3 (crosses), Xla-Xld (filled squares), 
and FS03L (plus sign). The lines show the best linear fits (Equation (26) for 
panel (b)). 

tion d23l in the second and third equalities, respectively. 

5.1.3. Outflow Flux 

We have seen in Section 14.11 and Figure |4{f) that MRI 
drives outgoing gas flow at the outer boundaries of the ac- 
tiv e layers. The MRI-driven outflow has been first observed 
bv lSuzuki & Inutsuka ( 2009) in shearing-box simul ations and 
been recently demonstrated by | Flock et all (1201 1 | ) in global 
simula tions. ISuzuki & Inutsukal (120091) and Suzuki et alJ 
( 12010b point out that this outflow might contribute to the 
dispersal of protoplanetary disks, although it is still unclear 
whether the outflow can really escape from the disks (see be- 
low). Meanwhile, MRI also contributes to the accretion of the 
gas in the radial direction. For consistent modeling of these 
two effects, we seek how the accretion stress and outflow flux 
are correlated with each other. 

We evaluate the o utflow mass flux in the following way. As 
seen in Section l4~Tl the temporally and horizontally averaged 
vertical mass flux (pv z ) is nearly constant at heights |z| > 5h. 
Using this fact, we define the outflow mass flux m w as the sum 
of | (pv z ) | averaged near upper and lower boundaries, 



X 



/'hiKl 



\(pvz)\dz 



(25) 



where /zb n d = 5\/2h « Ih is the height of the upper and lower 
boundaries of the simulation box. 

Figure fTOi a) shows m H , normalized by (p) m idC, versus a for 
all simulations. The di mensionless quantity m w /((p) m nc s ) is 
equivalent to C w used in lSuzuki et al.l d2010t) . For model Ideal, 
the value m „ = 5 x 10~ 5 (p) m j&c s is con sistent with the (3 z q = 10 5 
ideal run of Su zuki & Inut suka ( 2009). The dashed line shows 
the best linear fit m w /((p)^c s ) = 0.015a. It can be seen that 
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the linear fit captures a rough trend but still considerably over- 
estimates the outflow flux for models Ideal and Y4. As seen 
in Tabled these are the models in which a C0le dominates over 
otatm- This implies that the turbulence in the disk core (which 
is the source of a core ) does not contribute to the outflow. In 
Figure flTiT b). we replot the data by replacing a with a a tm- We 
find that m w more tightly correlates with a atm than with a. 
The best linear fit is found to be 



m w = 0.016a atm (p) m i d c J . 



(26) 



This result is consistent with the idea that the outflow 
is driven at the outer boundaries of the active layers 
(ISuzuki & Inutsu ka 2009), because the dominant contribution 
to a atm comes from heights very close to /jideai- 

Although outflow from the simulation box is a general phe- 
nomenon in our simulations, it is unclear whether the outflow 
leaves or returns to the disk. In fact, the outflow velocity ob- 
served in our simulations does not exceed the sound speed 
even at the vertical boundaries. Since the escape velocity 
is higher than the sound speed, this means that the outflow 
does not have an outward velocity enough to escape out of 
the disk. Acceleration of the outflow beyond the escape ve- 
locity has not b een directly demonstrated by previou s simu- 
latio ns as well (ISuzuki et al.ll2010t IFlock et alJl20lth . How- 
ever, ISuzuki et all (120101) point out a possibility that magne- 
tocentrifugal forces and/or stellar winds could accelerate the 
outflow to the escape velocity. If the escape of the outflow 
will be confirmed in the future, our scaling formula for m w 
will certainly become a useful tool to discuss the dispersal of 
protoplanetary disks. 

5.2. Saturation Predictors for the Accretion Stresses 

In the previous subsection, we have shown that the ampli- 
tudes of various turbulent quantities scale with the vertically 
integrated stresses a C ore and a atm . The next step is to find out 
how to predict a core and a atm in the saturated state from the 
vertical magnetic flux B z o (or equival ently /3-o) and the resis- 
tivity profile j]. As shown in Section 14721 the turbulent state 
of a disk depends on the resistivity only through the critical 
heights of the dead zone, h a and /i res . Furthermore, the val- 
ues of Iia and /i les are only weakly affected by the nonlinear 
evolution of MRI since the fluctu ation s in B z and p are small 
inside the dead zone (see Section |4~TT i. Therefore, we expect 
that the effect of the resistivity can be well predicted by the 
values of li\ and h Tes in the initial state, i.e., h\ t o and /i res ,o- 
With this expectation, we try to derive saturation predictors 
for a core and a atm as a function of f3 z o, /za.o, and /z les ,o- 

First, we focus on a core . Figure fTTT a) plots a core versus 
for all our simulations. We see that a core scales roughly 
linearly with /3Iq. The deviation from the linear scaling is 
expected to come from the difference in the dead zone size, 
i.e., /ia,o and /i res ,o- In Figure[l2] we plot the product a C oreA:0 
as a function of /z res .o- For models except Ideal and Y4, we 
find that atoreAfl is well predicted by a simple formula 



acore PzO = 510exp(-0.54/z res , //2). 



(27) 



Fi gure fTTT b) replot the data in Figure [TTJ a) by replacing P z q 

with exp(-0.54/i reSi o//i)/3ro- F° r models Ideal and Y4, Equa- 
tion (|27| | underestimates a core . As explained in Section l4~2l 
these models exhibit higher magnetic activity near the mid- 
plane than the other models because of no or a thin dead zone 
(2/za < h). We expect that the higher magnetic activity gives 
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Figure 11. Disk core accretion stress a C ore for all our simulations, (a) Vs. the 
inverse initial plasma beta fc 1 . (b) Vs. fc 1 multiplied by exp(-0.54/i lcs ,o//i) 
(see also Figure 1121 . (c) Vs. the final predictor function, Equation j28> 
(solid line). The symbols correspond to models Ideal (circle), X0-X3 (open 
squares), Y1-Y4 (triangles), W1-Y3 (crosses), Xla-Xld (filled squares), 
and FS03L (plus sign). The dashed lines in panels (a) and (b) are linear 
fits, shown only for reference. 




Figure 12. a CO re/3;0 versus /i res ,o for all runs. The symbols correspond to 
models Ideal (circle), X0-X3 (open squares), Y1-Y4 (triangles), W1-W3 
(crosses), Xla-Xld (filled squares), and FS03L (plus sign). T he s olid line 
shows an exponential fit for runs except Ideal and Y4 (Equation (27)). 
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Figure 13. Atmosphere accretion stress a dtm for all our simulations, (a) Vs. 
the inverse initial plasma beta /37 ' ■ (b) Vs. the final predictor function, Equa- 
tion )29t (solid line). The symbols correspond to models Ideal (circle), X0- 
X3 (open squares), Y1-Y4 (triangles), W1-W3 (crosses), Xla-Xld (filled 
squares), and FS03L (plus sign). The dashed line in panel (a) is a linear fit, 
only shown for reference. 

additional contribution to a core . Taking into account this ef- 
fect, we arrive at the final predictor function, 

a core = 510exp(-0.54/ 2res0 //i)/3 z "o +0.01 1 exp(-3.6/i A o/h). 

(28) 

Here, the numerical factors 0.011 and 3.6 appearing in the 
second term have been chosen to reproduce the results of runs 
Ideal and Y4, respectively. Figure [TTT c) compares the final 
fitting formula with the numerical data. It can be seen that 
Equation (1281 well predicts a core for all our models. Note that 
the second term of the predictor function is assumed to have 
no explicit linear dependence on f3~^ unlike the first term. In 
fact, it is possible to reproduce our data by multiplying the 
second term by a prefactor 3 x 10 s //3 z q. However, as we will 
see below, the absence of the prefactor makes the predictor 
function consistent with the results of ideal MHD simulations 
in the literature. 

The predictor function for a atm can be obtained in a similar 
way. Figure [T3f a) shows a atm versus /3" ' for all our runs. We 

find that a simple linear relation a atm = 530/3r ' well fits to the 
data except for models Ideal and Y4. This means that a aim 
is characterized only by /3" ' as long as the dead zone is thick 
(2h\ > li). To take into account the cases of thin dead zones, 
we add a term proportional to exp(-3.6/iA,o/^) as has been 
done for a C0le , and obtain 

a atm = 530^' +0.0043exp(-3.6/i A .o//i), (29) 

where the prefactor 0.0043 for the second term has been deter- 
mined to fit to the result of run Ideal. As seen in Figure [T3l b). 
Equation ( 1291 well predicts the value of a- dtm for all our runs. 

Our predictor functions indicate that the vertically inte- 
grated accretion stress is inversely proportional to /3-o when 
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Figure 14. Accretion stress (w.xy) versus gas density (p) at the upper bound- 
ary of the active layer (z = ^ideal) for all runs. The symbols correspond to 
models Ideal (circle), X0-X3 (open squares), Y1-Y4 (triangles), W1-Y3 
(crosses), Xla-Xld (filled squa res), and FS03L (plus sign). The solid line 
shows a linear fit (Equation (30}). 
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Figure 15. Vertically integrated accretion stress a versus the column density 
^active of the active layer for /3-o = 3 X 10 5 models. The symbols correspond 
to models Ideal (circle), X0-X3 (open squares), Y1-Y4 (triangles), W1-Y3 
(crosses), and FS03L (plus sign). The dashed line is a linear function, only 
shown for reference. 

a large dead zone is present (/ja > h). As we show below, this 
dependence originates from the magnitude of the accretion 
stress at the outer boundaries of the active layers, |z| w /zideai- 
When a dead zone exists, the dominant contribution to a 
comes from the accretion stress at that location (see Figures[5] 
and|6]). As shown Figure[l4j our simulations suggest that the 
accretion stress at |z| = /zideai obeys a simple relation 

(h^XAhmi) « 0. 18(p) (/lidealk 2 . (30) 

This means that the averaged accretion stress at |z| = /lideat is 
18% of the averaged gas pressure (pc 2 s ) at the same height. By 
the definition of /iid ea i, the gas density at |z| = /lideai is related to 
(B 2 ) at the same height as (p)(/zideat) = (27r) 2 (fi 2 )(/jid e af)/47rc 2 . 
Since our simulations suggest (-S 2 )(/j;deai) ~ 10-S 2 , the rela- 
tion means (p)(/!id ea i) ~ 10(27r) 2 B 2 /47rc 2 ~ 10 3 /3^'po- Using 
this fact, Equation (f30b can be rewritten into the linear relation 
between (wxy)(/zid e af) and p^: 

(w,v)(/! id eai)~100/3: Voc 2 . (31) 

When a dead zone is present, the level of a is determined by 
(w xv )(/z;deat) (see above), so we have a oc /37q. 

We remark that the vertically integrated stress does not 
scale linearly with the column density of active layers. This 
is shown in Figure Q3J where we compare a with the column 
density S a ctive of the active region h\ < \z\ < /zideai- We see 
that a decreases much more slowly than S a ctive when S a ctive is 
less than 10% of the total gas surface density. This reflects the 
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fact that the dominant contribution to a comes from the outer 
boundaries of the active zones, |z| w /i;deai- 

It is useful to see how the predictor functions work when 
a dead zone is absent. If h\$ = /z res .o = 0, Equations d28T i and 
d29l predict the total accretion stress a = 1 .0 x 10 3 /3~q +0.015. 
This implies that a is constant {a ps 10~ 2 ) for f3 z o > 10 5 and 

increases linearly with (a « 10 _2 (10 5 /Ao)) for /3 z0 < 10 5 . 
Strikingly, this pre diction is consistent with the finding by 
iSuzuki et al.1 (I20101 see their Figure 2). The existence of the 
floor value a w 10~ 2 at low net vertical magnetic fluxes (i.e., 
at high f3 z o) is also supporte d by recent stratifi ed MHD sim- 
ulations with zero net flux dDavis et alj 12010b . These facts 
suggest that our predictor functions are applicable even when 
a dead zone is absent. 

6. VERTICAL DIFFUSION COEFFICIENT 

As seen in Section |4] sound waves excited in the upper 
layers create fluctuations in the gas velocity near the mid- 
plane. It has been well known that fully developed MRI- 
driven turbulence causes the diffusion of sm all dust particles 
dJohansen & Klahil2005t iTurneret al.ll2010t) . However, it has 
not been fully understood how the sound waves propagating 
inside a dead zone affect the dynamics o f dust particles there. 
For example, ISuzuki & Inutsukal (12009b speculated that the 
sound waves might promote dust sedimentation by transfer- 
ring t he downward mome ntum to dust particles. On the other 
hand, iTurner et al.1 (120101) reported that the waves excite ver- 
tical oscillation of dust particles deep inside the dead zone 
and thus prevent the formation of a thin dust layer. Since dust 
sedimentation is crucial to planetesimal formation via gravita- 
tional instability, it is worth addressing here how it is affected 
by the velocity dispersion created by sound waves. 

Here, we focus on the dynamics of small dust particles, and 
model the swarm of the particles as a pass ive s calar as was 
previo usly done bv lJohansen & Klahrl d2005b and lTurner et al.l 
d2OT0l) . We assume that dust particles are so small and their 
stopping time r s is much shorter than the turnover time of 
turbulence (~ fT 1 ). We also assume that the dust density is 
lower than the gas density and hence the dust has no effect 
on the gas motion. Under these assumptions, the velocity of 
dust particles relative to the gas can be approximated by the 
terminal velocity Vr = -il 2 r s zi, where z is the unit vector for 
the z-direction. Then, the equation of continuity for dust is 
given by 

^+V-[p d (v+V r )] = 0, (32) 
at 

where pd is the dust density. Equation d32l has an advantage 
that the time step can be taken longer than t s in numerical 
calculation. 

We have solved Equation ( 1321 for four models (Ideal, XI, 
X3, and Yl) with the initial condition that the dust-to-gas 
mass ratio / = pd/p is constant throughout the simulation 
box. To extract the effect of the quasi-stationary turbulence, 
we insert the dust 100 (for models XI, X3, and Yl) or 250 (for 
model Ideal) orbits after the beginning of the MHD calcula- 
tions. The stopping time t s is set to t s = 0.01S1 -1 for model 
Ideal and t s = 0.001 SI -1 for the other models. The longer t s 
has been adopted for model Ideal to allow the dust to settle 
appreciably in the stronger turbulence. In reality, the stopping 
time of a dust particle depends on the gas density and hence 
on z, but we ignore this dependency for simplicity. 

Figure [16] shows the temporal evolution of the dust density 



at the midplane, pd.miA, for the four MHD runs. The solid 
curves show the horizontally averaged pd,m\A observed in the 
MHD runs. For comparison, the evolution of pd,miA in a hy- 
drostatic, laminar disk is also shown by the dotted curves. Ir- 
respectively of the presence or absence of a dead zone, the 
dust density observed in the MHD runs is higher than that in 
a laminar disk at all moments. This means that sound waves 
propagating in a dead zone do not promote but prevent dust 
settling as turbulence does in an active zone. 

To illustrate more clearly the diffusive nature of the ve- 
locity dispersion in a dead zone, we try to compare the 
above results with a simple advection-diffusion theory. Here, 
we consider a one-di mensional advection-diffusion equation 
dDubrulle et alj[l995b 



dpd _ d(p d V T ) | d f d_pd_ 
dt dz dz \ z 9z p 



(33) 



where Vt is the terminal velocity given above and D- is the 
vertical diffusion coefficient for dust. If the dust particles are 
sufficiently small (r s £l <C 1), D z is equal to the diffusion co- 
efficie nt for gaseous contaminants (e.g.. lYoudin & Lithwickl 
12007b . The first term in the right-hand side of Equation fl33T > 
represents the downward advection of dust due to settling, 
while the second term represents the diffusion of dust in 
the stratified gas. For disks with no or a small (\z\ < h) 
dead zone, it is known that Equation ( f33l > well describes the 
evolution of Pd if the diffusio n coefficient is assumed to be 
dFromang & Papa loizou 2006b 



D z ^{5v\)/n. 



(34) 



We here examine whether Equations (1331 and ( 1341 work 
well even when a dead zone is present. We solve Equa- 
tion d33l with D z = b(5v 2 )/fl, where the vertical distribu- 
tion of (5\> 2 ) is taken from temporally and horizontally av- 
eraged MHD data and b is a dimensionless fitting parameter. 
The dashed curves in Figure [16] show the predictions by the 
advection-diffusion model, where b is set to be 1.0, 0.5, 0.9, 
and 1.0 for runs Ideal, XI, X3, and Yl, respectively. It can 
be seen that the advection-diffusion model with b ~ 1 suc- 
cessfully reproduces the long-term evolution of the observed 
Pd.mid for all the models. It is striking that a constant b well re- 
produces the evolution of the dust density at all heights, as is 
shown in Figure[l7] In this figure, the solid and dashed curves 
show the vertical distribution of the dust-to-gas mass ratio 
f = pd/Pg observed in run XI and predicted by the advection- 
diffusion model with b = 0.5, respectively. In this run, the 
boundaries between the active and dead zones are located at 
|z| = /ja ~ 2.5/i (see Table [2J. However, Equation (l33l suc- 
cessfully predicts the evolution of / even if we do not change 
the value of b across the boundaries. This fact supports the 
idea that sound waves propagating across a dead zone con- 
tribute to the diffusion of dust particles just as turbulence does 
in active zones. 

It is worth mentioning here that the diffusio n coefficient 
D 7 inc rease s with \z\ as has been poin ted out by ITurner et al.l 
d2006b and iFromang & Nelson! (I2009b . This effect is partic- 
ularly significant at high altitudes where the gas density is 
much lower than that at the midplane, because D z oc (5v 2 ) is 
roughly proportional to the inverse of the gas densityQ The 
dotted curves in Figure (T% show how Equation d33l l would 



4 Fromane & Nelson (2009) used a diffusion coefficient quadratic in z to 
explain the dust distribution in their MHD simulation. Our finding D z oc p 
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Figure 16. Temporal evolution of the dust density at the midplane in various MHD runs. The four panels show the results for models Ideal (upper left), XI (upper 
right), X3 (lower left), and Yl (lower right). The solid curves show the horizontally averaged dust density pd.mii observed in the MHD runs, while the do tted curve 
show the evolution of pj m ii in a laminar disk. The dashed curves are the predictions from the one-dimensional advection-diffusion equation (Equation j33t ). The 
diffusion coefficient D- in Equation j33t is assumed to be proportional to (<5v?)/n, and the proportionality factor has been chosen to best reproduce the evolution 
of Prf m id in the MHD runs. 
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Figure 17. Snapshots of the vertical distribution of the dust-to-gas mass ra- 
tio /= p d /p at t = 100, 250, 400, and 550 orbits for model XI. The solid 
curves show the horizontally averaged / observed in the MHD simulation. 
The dashed curves show the sol utions to the one-dimensional advection- 
diffusion equation (Equation (33)) with D.(z) = 0.5(5vJ)(z)/n. The dotted 
curves show the solution to Equation j33> with a constant diffusion coe£5- 
cientZ) ; = 0.5{<5v?) mld /£l 

fail to predict dust evolution if one assumed a constant dif- 
fusion coefficient D z = 0.5(<5v?) m jd/ll We see that the con- 
stant diffusion coefficient model significantly underestimates 
the dust density at |z| ^> h. This fact will merit considera- 

does not contradict their assumption because p~' oc exp(z 2 /2h 2 ) m 1 +z 2 /2h 2 
near the midplane. 



tion when modeling the chemical evolution of protoplanetary 
disks, in which the vertical m ixing of molecules is of impor- 
tance dHeinzeller et al.ll201 lb . 

Finally, we give a simple analytic recipe for the vertical dis- 
tribution of D z . It is useful to rewrite Equation J34t in terms 
of (5v 2 ), for which the scaling relation (Equation d22l i) and 
predictor function (Equation d28l l) are available. Table 2 lists 
the ratio of (<5v?) m id to (<5v 2 ) m jd for all our simulations. It can 
be seen that (<5v 2 )mid « (32% ± 15%) x (<5v 2 ) m id, indicating 
that (5v 2 ) m id is roughly equal to a third of ((5v 2 ) m ;d. Fur- 
thermore, the ratio (i5v 2 )/(<5v 2 ) is approximately constant in 
the disk core, as is illustrated in Figure |4jc). Based on these 
facts, we approximate (<5v 2 ) as (Sv 2 ) /3 in the disk core. Using 
this approximation together with the scaling relation for (<5v 2 ) 
(Equation d22l). we rewrite Equation ( 1341 as 

«0.3(a coreC 2 /ft)exp^|^. (35) 

If one uses this equation together with the predictor function 
for a cor e (Equation d28l)). one can calculate the vertical distri- 
bution of D z in the disk core for given /3 z q and r\. 

7. DISCUSSION: EFFECTS OF NUMERICAL 
RESOLUTION 

All MHD simulations presented in the previous sections 
were performed with the numerical resolution of 40 x 80 x 
200 grid cells for the simulation box of size 2V2h x 8V2/1 x 
10V2/L Here, we examine how the numerical resolution af- 
fects the saturated state of turbulence. 
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Figure 18. Saturated values of various quantities versus numerical res- 
olution for model XI. From top to bottom: a, a CO re, (5v 2 }mid/ c ?i 
(<5p 2 )mid/(p)mid' an d ™w / (p) mid£s ■ The horizontal axis shows the number 
of grid cells per length V2h in the vertical direction, n z . The value n z = 20 
corresponds to the resolution adopted in this study. 

We carry out XI simulations with changing the numerical 
resolution to 20 x 40 x 100 cells and 80 x 160 x 400 cells. 
Figure Q~8] compares the saturated values of various quantities 
obtained from the two runs with the values from the original 
XI run (Table 2). Here, the horizontal axis shows the number 
of grid cells per length y/2h in the vertical direction, n z (see 
footnote [TJ. The value n z = 20 corresponds to our original 
resolution. We see that the change of the resolution hardly af- 
fects the integrated accretion stresses a and a core and outflow 
flux m w , suggesting that the resolution of n z = 20 is sufficient 
for these quantities to converge well. By contrast, the veloc- 
ity and density dispersions (Sv 2 ) m id and (Sp 2 ) m id increase with 
improving the resolution. Since the energy input rate to tur- 
bulence should be the same if the integrated accretion stress 
is unchanged, the resolution dependence of the velocity and 
density dispersions is expected to mainly come from artificial 
dissipation of sound waves in the simulation box. However, 
we also see that this effect becomes less significant as the res- 
olution is improved. Detailed inspection shows that the frac- 
tional increase in (<5v 2 ) m jd is 81% when going from n z = 10 to 
n z = 20 but is 40% when going from n z = 20 to n z = 40. This 
suggests that the amplitudes of the velocity and density fluc- 
tuations should converge to finite values in the limit of high 
resolutions (n z — > oo). This is to be expected, since sound 
waves in a stratified disk physically dissipate through, e.g., 
shock formation, particularly at high altitudes where the gas 
density is low and therefore the amplitudes of the waves be- 
come large. We find that the data for (Sv 2 ) m i d /c 2 shown in 
Figure QJ lie on a curve (Sv 2 ) mid /c 2 = 0.012- 0.0015n7 015 . 
This implies a converged value of (Sv 2 ) m a/c 2 « 0.012, which 
is five times higher than that obtained in our n z = 20 simula- 
tion ((5v 2 )mid/Cj « 0.0024). From this estimate, we see that 
(<5v 2 ) and (Sp 2 ) could be underestimated by a factor of several 
in the simulations presented in this study. 

In summary, we find that the outflow mass flux and verti- 
cally integrated accretion stress converge well within our nu- 
merical resolution. This suggests that the predictor functions 
for acore and a- dtm (Equations (l28l and d29| i) and the scaling 
relation between m w and a atm are hardly affected by the res- 
olution. On the other hand, the amplitudes of sound waves 
could be underestimated by a factor of several because of the 
finite grid size. Future high-resolution simulations will enable 
to better quantify the scaling relations between (i5v 2 ) and a core 
and between (Sp 2 ) and ow (Equations (l22l and (T24l). 



8. SUMMARY 

Good knowledge about the turbulent structure of protoplan- 
etary disks is essential for understanding planet formation. 
To provide an empirical basis for modeling the coevolution 
of dust and MRI, we have performed MHD simulations of a 
vertically stratified shearing box with an MRI-inactive "dead 
zone" of various sizes and with a vertical magnetic flux of 
various strengths. Our findings are summarized as follows. 

1. We have introduced the critical heights (h ic \ e . d \, h\, and 
/ires) that characterize the MRI in a stratified disk (Sec- 
tion |3). We have found that the vertical structure of 
MRI-driven turbulence depends on the resistivity pro- 
file only through the critical heights for the dead zone 
(fiA and hies) and is i nsens itive to the detail of the resis- 
tivity profile (Section l4T2l . 

2. In the "disk core" (|z| < ZiideaiX the density-weighted ve- 
locity dispersion (p) (8v 2 ) is ne arly constant along the 
vertical direction (Section 14. 11 1. This means that the 
velocity dispersion is approximately inversely propor- 
tional to the gas density. Weak dependence on z is also 
found for (Sp 2 ) / (p), meaning that the density fluctua- 
tion (Sp 2 ) 1 / 2 is proportional to the square root of the 
averaged density. 

3. The accretion stresses in the disk core and "atmo- 
sphere" (\z\ > /iideai) differently co ntrib ute to the tur- 
bulent structure of a disk (Section |5. 11 1. The velocity 
dispersion (Sv 2 ) and density fluctuation (Sp 2 ) in the 
disk core depend linearly on the accretion stress inte- 
grated over the core, a core (Equations (fSTJ and d23l). By 
contrast, the outflow mass flux m u . depends linearly on 
the stress integrated over the atmosphere, a a tm (Equa- 
tion (1261). 

4. We have obtained simple empirical formulae that pre- 
dict the vertically integrated stre sses a core and a at m in 
the saturated state (Section 15.21 Equations (l28l and 
(|29]>). These are written as a function of the strength 
of the vertical magnetic flux (or /3-o) and the critical 
heights of the dead zone measured in the nonturbulent 
state (/i res ,o and /ia,o). These predictor functions to- 
gether with the saturation relations described above al- 
low to calculate various turbulent quantities for a given 
resistivity profile and a net vertical flux. 

5. We have confirmed that the vertical diffusion coeffi- 
cient D z of contaminants is given by D z w (Sv 2 ) /Q both 
inside and outside a dead zone (Section|6]l. This implies 
that sound waves propagating across a dead zone con- 
tribute to the diffusion of dust particles just as turbu- 
lence does in active zones. We have obtained a simple 
analytic recipe for the vertical distribution of D z as a 
function of a core on the basis of our MHD simulation 
data (Equation (l35l>). 

The empirical formulae obtained in this study enable us to 
predict the amplitudes of various turbulent quantities in a pro- 
toplanetary disk with a dead zone. The steps to be performed 
are as follows. 

1. Prepare the vertical profile of the ohmic resistivity rj, 
and find h\ and /i les from Equations ( TTZb and (fT3l . 
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A realistic profile of rj in the presence of dust parti- 
cles can be obtained by solving the ionization state of 
the g as and the charge state of the dust simul taneously 
(e.g. ISano et al.l2000llllgner & Nelsonl2006tlOkuzumil 
120091) . 

2. Calculate a core and a atm using the predictor functions, 
Equations @ and d2% 

3. One can now calculate the turbulent viscosity of 
the disk as = (3 /2)(q core + a a tm)c 2 / O (see Equa- 
tions ( TT~8T >. (fT~9T >, and (l20ii). The vertical distribution of 
the gas velocity dispersion and density fluctuation in 
the disk core (|z| < /ijdeai) can be calculated from Equa- 
tions d22l i and d24b . respectively. The outflow mass flux 
can be evaluated from Equation ( f26b . For the diffusion 
coefficient in the disk core, one can use Equation ( [35l . 

When using our empirical formulae, it should be kept in 
mind that our scaling relations for velocity and density fluc- 
tuations (Equation (l22l i and d24l i) could underestimate their 
mean-squared amplitudes relative to the integrated accretion 
stress by a factor of several because of the numerical dissipa- 
tion of sound waves (Section|7]i. Future high-resolution sim- 
ulations will allow to better quantify the saturation level of 
sound wave amplitudes. 

We are grateful to Neal Turner and Takayoshi Sano for pro- 
viding us with their MHD simulation data that have motivated 
us to start this study. We also thank Shu-ichiro Inutsuka, 
Takeru Suzuki, Taku Takeuchi, Hidekazu Tanaka, Takayuki 
Tanigawa, Mordecai-Mark Mac Low, and the anonymous ref- 
eree for useful discussion and fruitful comments. Calculations 
were made on the Cray XT4 at the CfCA, National Astronom- 
ical Observatory of Japan, and SRI 6000 at the Yukawa Insti- 
tute for Theoretical Physics, Kyoto University. S.O. is sup- 
ported by a Grant-in- Aid for JSPS Fellows (22 ■ 7006) from 
the MEXT of Japan. 



in the gas phase, «, = n e . Equation ( lAlb leads to the electron 
abundance 

x e = — = a — ^— oc VCexp ( -^—] , (A2) 

which means that the resistivity is proportional to 
£~ l//2 exp(-z 2 /4/i 2 ). Thus, if the vertical dependence of 
( can be neglected, the resistivity profile is given by Equa- 
tion © with h n = y/2h. Note that Equation © is derived 
instead of Equation dU if cosmic-ray ionization is assumed 
and the attenuatio n of cosmic rays towar d the midplane is 
taken into account dFleming & Stonell2003l) . 

If the total surface area of dust particles is large, recombina- 
tion occurs mainly on dust surfaces. In this case, the equation 
for the ionization-recombination equilibrium is given by 

C"« = Idendtle, (A3) 

where jde is the sticking rate coefficient for dust-electron col- 
lision and rid is the number density of dust particles. The 
sticking rate coefficient depends on the charge of the dust 
particles, and u ltimately on n e via the charge neutrality (see 
IOkuzumi| [2009). but we will i gnore this dependence in the 
following. From Equation ( IA3b , we have 

x e = — — oc Cexp (^jj) , (A4) 

where we have assumed that rid oc exp(-z 2 //i^) with hd be- 
ing the scale height of the dust particles (in fact, one can 
show using Equation d33l that rid obeys a Gaussian distribu- 
tion in sedimentation-diffusion equilibrium if D z oc t s oc n~ l oc 
exp(z 2 /2/z 2 ); the condition t s oc n~ l is satisfied if the size of 
the dust particles is smaller than the mean free path of the 
gas). Thus, ignoring the dependence of ( on z, the resistivity 
rj oc x~ l is given by Equation (HJi with h v = hd. 
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APPENDIX 

IONIZATION DEGREE AND OHMIC RESISTIVITY 
IN PROTOPLANETARY DISKS 

In this Appendix, we explain how the resistivity profile 
adopted in this study (Equation (@]i) is related to realistic re- 
sistivity profiles in protoplanetary disks. Since the resistivity 
is inversely proportional to the io nization degree (more pre- 
cisely, the electron abundance; see lBlaes & Balbuslll994l) . we 
will see how the ionization degree depends on the height z 
above the midplane. 

Recombination occurs in the gas phase and on dust sur- 
faces. The gas-phase recombination dominates if the total 
surface area of dust particles is negligibly small. In this case, 
the equation for the ionization-recombination equilibrium is 
given by 

(n„ = iieriin e = j ie n 2 e , (Al) 

where £ is the ionization rate (the probability per unit time at 
which a molecule is ionized), «„, «, , and n e are the number 
densities of neutrals, ions, and electrons, respectively, and 7^ 
is the gas-phase recombination rate coefficient. The second 
equality in the above equation assumes the charge neutrality 
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